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Summary. — During the last 15 years, there has been much progress in defining the 
nuclear Hamiltonian and applying quantum Monte Carlo methods to the calculation 
of light nuclei. I describe both aspects of this work and some recent results. 



1. Introduction 

The goal of ab-initio light-nuclei calculations is to understand nuclei as collections 
of nucleons interacting with realistic (bare) potentials through reliable solutions of the 
many-nucleon Schrodinger equation. Such calculations can study binding energies, exci- 
tation spectra, relative stability densities, transition amplitudes, cluster-cluster overlaps, 
low-energy astrophysical reactions, and other aspects of nuclei. Such calculations are also 
essential to claims of sub-nucleonic effects, such as medium modifications of the nuclear 
force or nucleon form factors; if a reliable pure nucleonic degrees of freedom calculation 
can reproduce experiment then there is no basis for claims of seeing sub-nucleonic de- 
grees of freedom in that experiment (beyond the obvious fact that the free-space nucleon 
interactions are a result of sub-nucleonic degrees of freedom). 

There are two problems in microscopic few- and many-nucleon calculations: 1) de- 
termining the Hamiltonian, and 2) given H, accurately solving the Schrodinger equation 
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for A nucleons; I will discuss both of these in this contribution. The two-nucleon (NN) 
force is determined by fitting the large body of NN scattering data. Several modern NN 
potentials are in common use. The Argonne vig is a local potential written in opera- 
tor format; this potential is used in the calculations described here, and is presented in 
some detail below. Other modern potentials are generally non-local; some of them are 
discussed in other contributions to this school. 

It has long been known that calculations with just realistic NN potentials fail to re- 
produce the binding energies of nuclei; three-nucleon (NNN) potentials are also required. 
These arise naturally from an underlying meson-exchange picture of the nuclear forces 
or from chiral effective field theories. Unfortunately, much NNN scattering data is well 
reproduced by calculations using just NN forces, so the NNN force must determined 
from properties of light nuclei. In this contribution the recent Illinois models with 2tt 
and 3ir rings are used. 

Our understanding of nuclear forces has evolved over the last 70 years: 

• Meson-exchange theory of Yukawa (1935) 

• Fujita-Miyazawa three-nucleon potential (1955) 

• First phase-shift analysis of NN scattering data (1957) 

• Gammel- Thaler, Hamda- Johnston and Rcid phenomenological potentials (1957- 
1968) 

• Bonn, Nijmcgcn and Paris field-theoretic models (1970s) 

• Tuscon-Melbourne and Urbana NNN potential models (late 70's-early 80's) 

• Nijmegen partial wave analysis (PWA93) with x 2 /dof^ 1 (1993) 

• Nijm I, Nijm II, Reid93, Argonne vi 8 and CD-Bonn (1990s) 

• Effective field theory at N 3 LO (2004) 

References for a number of these developments are given in the following sections. 

Accurate solutions of the many-nucleon Schrodinger equation have also evolved over 
many decades: 

• 2 H by numerical integration (1952) - a pair of coupled second-order differential 
equations in 1 variable. At the time this took "between 5 and 20 minutes for the 
calculation and the printout another 5 minutes" [T]! 

• 3 H by Faddeev (1975-1985) 

• 4 He by Green's function Monte Carlo (GFMC) (1988) 

• A = 6 by GFMC and No-core shell model (NCSM) (1994-95) 

• A = 7 by GFMC and NCSM (1997-98) 

• A = 8 by GFMC and NCSM (2000) 

• 4 He benchmark by 7 methods to 0.1% (2001) 

• A = 9, 10 by GFMC and NCSM (2002) 

• 12 C by GFMC and NCSM (2004-) 

• 16 O by Coupled Cluster (CC) (2005-) 
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References for the A=3,4 calculations may be found in Ref. [2]; the GFMC calculations 
are the subject of this paper; the NCSM are discussed in Petr NavratiPs contribution to 
this Course; and CC results may be found in Ref. [3]. 

This contribution is limited to Variational Monte Carlo (VMC) and GFMC calcula- 
tions of light nuclei. Section [2] describes the Hamiltonians used and sections [3] through [5] 
describe the computation methods. Section [6] gives a number of results for energies of 
nuclear states; Sec. [7] describes GFMC calculations of scattering states; and Sec. [8] gives 
some results for densities. Finally some conclusions and prospects for the future are 
presented in Sec. [9j 



2. Hamiltonians 

The nuclear Hamiltonian used here has the form 

(1) a »'< • X>j • E v ^ ■ 

i i<j i<j<k 

Here Ki is the non-relativistic kinetic energy, including m n — m p effects, Vij is the 
AW potential and Vijk is the NNN potential. 

2T. Argonne v^. - The AW potential (vij) is Argonne vl8 [3] (AV18) which has the 
form 

(2) ^=4+^+^ + ^ 7S . 

The vjj is a very complete representation of the pp, pn and nn electromagnetic terms, 
including first- and second-order Coulomb, magnetic, vacuum polarization, etc., compo- 
nents with form factors. (Ref. [5] provides a heuristic introduction to AV18.) 

The vfj is the isoscalar one-pion exchange potential represented as a local operator: 



('-!) <V, \ %^ ' X, r,-r, . 



(4) Xy = T(m n rij) Sij + Y{m^r i:j ) ai ■ <Tj , 

(5) Y(x) = — £(r) , 

x 

(6) T(x) = + I + l) Y(x) £(r) , 

\x z x ) 

(7) eM = (l-e^ 2 ). 



where r^, <7j and Sij are isospin, spin and tensor operators, respectively. In light nuclei, 
(vfj) contributes ~85% of (vy). 

The remaining isospin-conserving terms are 

(8) 4 = £ v virij)0% , 

p=l,14 
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Fig. 1. - Two-pion exchange terms in the Illinois NNN potentials. 



(9) Of" 1 ' 14 = [1, a z ■ a j , Sij , L ■ S, L 2 , L 2 ^ ■ ^, (L • S) 2 ] ® [1, n ■ Tj ] , 

where v p (r) has short-, intermediate-, and long-range components. The long-range com- 
ponents are just the Y(r) and T(r) of the one-pion potential and are present only for 
those operators that have contributions from one-pion exchange. The intermediate-range 
components are proportional to T 2 (r) and the short-range component is of the Woods- 
Saxon form. 

Finally, v^ IB is the strong charge independence breaking part of the potential and 
consists of four operators: 

(10) O^ 15 ' 18 = [1, fa • a,), ® T tj , (r zl + r zj ) . 

The long-range part of O p=15,17 comes from one-pion exchange by inserting m„+- or 
m^o in Eqn. ([3]and[l]) and using f 2 NN oc m„. 

The parameters in the short- and intermediate-range components were determined by 
making a direct fit to the 1993 Nijmegen data base [6J [7] containing 1787 pp and 2514 np 
data in the range — 350 MeV, the nn scattering length, and deuteron binding energy. 
The fit of approximately 40 parameters results in a x 2 /d.o.f. of 1.09, which is typical of 
1990's NN potentials. 

2'2. Illinois Vijk- - The three-nucleon potential used for most of the examples pre- 
sented here is the Illinois-2 [8]. It consists of two- and three-pion terms and a simple 
phcnomcnological repulsive term: 

(11) = v$ + v$ + v$k ■ 

The two-pion term, illustrated in Fig.[TJ contains P- and 5*-wave 7r^V-scattering terms: 

(12) V 2 l = V 2 Z' P + V 2 l' S • 

V / tjk ijk 1 i]k 

The P-wave term (left panel of Fig. [1]) is the well-known Fujita-Miyazawa [9l [10] term 
which is present in all realistic NNN potentials. It has the form 

( 13 ) V 1]k F = A ^,P H \ X :r X ■ T J, T 3 ■ T k} + -^Xij,X jk \[Ti ■ Tj,Tj ■ T k ] , 

cyclic 
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Fig. 2. - Three-pion ring terms in the Illinois NNN potentials. 



where Xy is defined in Eq. (01 . This is the longest-ranged nuclear NNN potential and 
is attractive in all nuclei and nuclear matter. However it is very small or even slightly 
repulsive in pure neutron systems. 

The second panel of Fig. [T] represents the ,5- wave part of V?£. This term was intro- 
duced in the Tuscon-Melbourne NNN potential |TTj and is required by chiral perturbation 
theory. However, in practice it is only 3%-4% of V?£' in light nuclei. 

The three-pion term (Fig. \2§ was introduced in the Illinois potentials. It consists of 
the subset of three-pion rings that contain only one A mass in the energy denominators. 
Even so it has a quite complicated form which is given in Ref. [S]. An important aspect 
of this structure is that there is a significant attractive term which acts only in T—3 /2 
triples. In most light nuclei (V$.) < O.l(F^) 

The final term in the NNN potential, V^ k , represents all other diagrams including 
relativistic effects. It is strictly phenomenological and purely central and repulsive: 

(14) V* k = A R ]T T\m v r ij )T' z {m v r jk ) . 

cyclic 

This repulsive term is principally needed to make nuclear matter saturate at the proper 
density instead of a too-high density and to obtain a hard enough equation of state for 
neutron matter. 

The coupling constants ^^p, A^, and Ar were adjusted to fit 17 nuclear levels for 
A < 8. The V^£' is too weak to be determined by fitting and its coupling was left at 
the value predicted by chiral perturbation theory. 

In light nuclei we find 

(15) (V ijk ) ~ (0.02 to 0.09) (%) ~ (0.15 to 0.6) (if) 

where the large fraction of (H) is due to a large cancellation of K and i?y . From this we 
expect 

(16) (V 4N ) ~ 0.06(F ijfc ) ~(0.02to0.04)(ii)~(0.5to2.)MeV. 

This is comparable to the accuracy of our calculations. Even if more accurate calculations 
could be made, it would probably not be possible to disentangle four-nuclcon potential 
effects from uncertainties in the fitted parameters of Vijk- 
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Fig. 3. - GFMC computations of energies for the AV18 and AV18+IL2 Hamiltonians compared 
with experiment. 



2'3. What Makes Nuclear Structure? . - We have defined a very complicated nuclear 
Hamiltonian and it is reasonable to ask if it is all necessary to reproduce the structure of 
light nuclei. A study [T2] was made of this in which features of the nuclear Hamiltonian 
were systematically removed and the effects on nuclear level energies investigated. For 
each simplification of the two-nucleon part of H, the remaining terms were readjusted 
to continue reproducing as many low partial-wave phase shifts, and the deuteron, as 
possible. 

Figure [3] shows the energies of various nuclear states. For each isotope there are three 
sets of energies: the right-most are the experimental values, the left-most are the results 
of the GFMC calculations to be described using just the AW potential AV18, and the 
middle ones are GFMC calculations using the AV18+IL2 Hamiltonian. The AV18+IL2 
results are generally in good agreement with the data; the rms deviation is ~0.75 MeV. 
However without the IL2 NNN potential the comparison to data gets steadily worse as 
the number of nucleons increases. This is a general result that has also been obtained 
by others using different many-body methods and different NN potentials. 

Figure [4] shows the effects of making further simplifications to H beyond removing 



Quantum Monte Carlo Calculations of Light Nuclei 



7 




What Makes 2s[111 
Nuclear Level Structure? 
GFMC Calculations 



AVI' AV4' AV6' AV8' IL2 



7 Li 



L S[4], 
'D[4] 
'G[4]l 
3 D[31}= : 
3 P[31]# 
5 D[22]f 
1 S[22]f • 
3 P[211] 



8 B( 



Fig. 4. - Nuclear energy levels for various simplifications of the Hamiltonian 



the NNN potential. Here the right-most results are again for the full AV18+IL2 Hamil- 
tonian and thus are close to the data. The next set of results to left are for the AV8' 
AW potential [T3] with no NNN force. With eight operators ([1, cTj • (Xj, Sy,L • S] ® [1, r, ■ 
Tj]), AV8' can reproduce AV18 results for eight partial waves; these are chosen to be ^So, 
3 5i, 3 Di, ei, x Pi, and 3 Po,i.2> (ei is the 3 Si- 3 Di mixing angle). (Strictly speaking, the 
3 i-2 potential of AV18 is reproduced but, because the 3 i<2 and €2 potentials are different 
from those of AV18, the 3 P2 phase shifts are not reproduced.) This potential is more 
attractive in nuclei than AV18 and more than makes up for the lost binding due to the 
removal of the NNN potential. In general it gives a good qualitative picture of nuclear 
energies. 

To the left of the AV8' results are results for AV6' which does not have L-S terms. 
Not surprisingly, these have only negligible spin-orbit splittings. In addition 6,7 Li are 
essentially unbound to breakup into a+d or a+t (the dashed lines show the indicated 
thresholds for each Hamiltonian). The next simplification is AV4' which contains no 
tensor force. With this force, the deuteron no longer has a D state but still has the 
correct binding energy. The ISo, a £i, ^i, and an average 3 Pj partial waves of AVI 8 are 
reproduced. This simplified force results in spurious degeneracies of nuclear levels and 
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somewhat overbinds all the shown nuclei. In particular 8 Be is bound against breakup 
into two alpha particles - an important failure because bound 8 Be would result in very 
different stellar nucleosynthesis. 

Finally the left-most results are for AVI', a pure central force that is an average of 
the ^So and 3 Si potentials of AV18. This produces an "upside-down" spectrum in which 
states with the least spatial symmetry are most bound. More importantly, there is no 
nuclear saturation; each increase in A results in much more binding and there are no 
A=5,8 mass gaps which are essential to big-bang nucleosynthesis and stellar evolution. 

Ref. [12] shows results for several other Hamiltonians including an AV2' that contains 
only central and space-exchange terms and thus is very similar to the popular Volkov 
potentials [2]. Besides erroneously binding the dineutron, this potential has the strange 
feature of binding 6 He but not 6 Li so that A=Q beta decay would be in the wrong 
direction. The conclusion of this study is that one needs almost the full, complicated, 
Hamiltonian to do realistic nuclear physics. 

3. Quantum Monte Carlo Methods 

The many-body problem with the full Hamiltonian described above is very difficult 
as is indicated by the slow progress over the last half-century that is outlined in the 
introduction. We need to solve 

(17) H ^(n,^,- ■ ■ ,?a;si,S2,- ■ ■ ,SA;h,t2,- ■ ■ ,tA) 



equations in 3 A — 3 variables (the number of isospin states can be reduced; see Sec l4'2[) . 
For 12 C this corresponds to 270,336 coupled equations in 33 variables. 

Furthermore, the coupling is strong; the expectation value of the tensor component 
of v* [Eq. [3] is approximately 60% of the total (uy ) but it is identically zero if there are 
no tensor correlations. Thus we cannot perturbatively introduce the couplings. 

We use two successive quantum Monte Carlo (QMC) methods to solve this problem. 
The first is variational Monte Carlo (VMC) in which a trial wave function, containing 
variational parameters, is posited and the expectation value of the Hamiltonian com- 
puted using Monte Carlo integration. In practice we have not been able to formulate 
accurate enough trial wave functions and so the second step, Green's function Monte 
Carlo (GFMC), is needed to iteratively project the exact eigenfunction out of the trial 
wave function. These two methods are described in the following sections. School and 
review articles on these methods are Refs. [HJ HH \T7\. Detailed descriptions may be 



E ^(fx,^, - ■ ■ ,r , A;s 1 ,s 2 r ■ ■ ,s A ;t 1} t 2 r ■ ■ ,tA) , 




found in Refs. QH EH OH CM [21] • 
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4. Variational Monte Carlo 

In VMC we start with a trial wave function, ^y, which contains a number of vari- 
ational parameters. We vary these parameters to minimize the expectation value of 
H, 

(18) (* r |* r ) ~ E °- 

As indicated, the resulting Et is, by the Raleigh-Ritz variational principle, greater than 
the true ground-state energy for the quantum numbers (J 7r , J z , T, and T z ) of 'J't- A 
simplified form of our trial wave functions is 

(19) |tf T > = + U l3 + Z k U tjk )]Y[f c (n 3 )\<I>) . 

Here / c (r) is a central (mostly short-ranged repulsion) correlation, Uij are non-commuting 
two-body correlations induced by Vij, and Uijk is a simplified three-body correlation from 
V ijh . 

More specifically, 



(20) = u p(nj)0: 



v 

ij 

p=2,6 



contains Tj • Tj, Oi ■ o~j, oi ■ o~j Ti ■ Tj, Sij, and Sy Tj • r operators, of which the Sy Tj • r is 
most important due to the already noted strong tensor contribution from v* . The / c (r) 
and u p (r) are solutions of coupled differential equations with as input [15] . 

The $ (see below) is fully antisymmetric; hence the rest of Eq. (fH)|) must be symmet- 
ric. But the Uij do not commute; for example 

(21) [<Ji ■ <7 2 , <7i ■ er 3 ] = 2i <7i ■ (o- 2 x er 3 ) . 

The symmetrizer S fixes this by summing over all [ ^^l -1 ^ ]! permutations of the ordering 
in rii<j- I n practice this is done by using just one Monte Carlo chosen ordering per wave 
function evaluation. 

4T. The one-body part of^r, _ The one-body part of \&t, is a lfuj shell-model 
wave function. It determines the quantum numbers of the state being computed and is 
fully antisymmetric. For 3 H and 3 ' He, $ can be antisymmetrized in just spin-isospin 
space, for example 

(22) m 3 H,Mj = !)) = ^(|pT nt ni) - bT raj n]) + |n|pT »T) 

- |ratPt raj.) + |raf raj p|) - |ra| rafpT) ) ■ 



10 



Steven C. Pieper 



For A > 4 we need P-wave radial wave functions in order to antisymmetrize <E>; the 
antisymmetrization is achieved by summing over all partitions of the A nucleons into 
four S'-shell nucleons (the a core) and A-4 P-shell nucleons which are antisymmetrically 
coupled to J 7 ' and T. To make $ translationally invariant, we express all functions of 
single-particle positions as functions of position relative to the center of mass of the A 
nucleons or of some sub-cluster of them. The one-body wave functions are solutions 
of Woods Saxon potentials containing several variational parameters. If desired the 
separation energy of the these one-body wave functions can be fixed at the experimental 
value to guarantee that the ^ T has the correct asymptotic form. In general <& has 
several spatial-symmetry components depending on how many ways a state of the desired 
quantum numbers can be constructed in the P-shell basis. For example the 6 Li $ have 
the form 



(23) |$) = AY,Pls\MLSJMTT 3 ) 1234 .. 56 ) 



LS 



(24) <S> e (LSJMTT 3 ) 123i .. 56 = $ 4 (OOOO) 123 40p i (P tt5 )^p (Re*) 

^[Y lmi (tt a5 )Y lm ,(n a6 )] LML x [X5(^m s )x6(^m' s )] S M s 

x M^ 3 K(^ 3 )]tt 3 • 



JM 



A Ihw LS'-basis diagonalization determines the 0ls- 



4'2. Representing in the computer. - The wave function, ^t^i, ^2, • • • , ?a), is a 
complex vector in spin-isospin space with dimension [N$ components for spin] x [Nt com- 
ponents for isospin]. The number of spin states is 2 A . However for even A, if we choose 
to use Mj = 0, we can calculate and retain only half the spin vector, say that part with 
positive spin for the last nucleon, and obtain the other half of the vector by time-reversal 
symmetry. The number of isospin states, Nt, depends on the isospin basis being used: 



(25) N T = for a proton-neutron basis , 

(26) = A/ ,2 +~r + i ( A/2 + t) for a good isos P in basis • 



Potentials (vij, Vijk) and correlations (tiy, f/yfc) involve repeated operations on \&. 
For example Oi ■ aj may be written as 



(27) 
(28) 



Oi ■ o~j = 2(cr+o- j + a l of) + o\o] 
1 



2P" 
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(29) 
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Here Pfj exchanges the spin of i and j. Consider the spin part of an A=3 wave function; 
<Tj • (ij will not mix different isospin components and, for different i and j, will separately 
act on different, non-contiguous, 4-element blocks of 



(30) 



* = 



(31) 



cr 2 • a 3 ^/ = 



( a TTT \ 
a TTi 
a TIT 
a tU 
a lTT 
a lTl 
a UT 

V a ni ) 



HT1 



2a 



ur 



l TTl 



2am - a TlT 

a lTT 
2a UT - a lTl 
2a lTl _ a UT 



a i • 0"2 V I / = 



0"3 • = 



a TTl 
2a m - a T | T 

2a lti ^ a TH 
2 «Tlt _ a lTT 

2a ui _ a m 



2a 



^TTT 
HTT ^ a TU 

2a UT _ a TH 
2a m - a iTT 

2a TU _ a llT 



Similarly, the tensor operator is 



(32) Sij = 3 (Tj • fij (Tj 



(33) 



/ z 2 - 1/3 
z(x + iy) 
z(xiy) 

V ( x + w) 2 



z(x - iy) 
-z 2 ~l/3 
x 2 + y 2 -2/3 
—z(x + iy) 



z(x - iy) 
x 2 + y 2 -2/3 
-z 2 + l/3 
—z(x + iy) 



{x-iyf \ 
-z(x - iy) 
~z(x - iy) 
z 2 - 1/3 / 



where etc. As shown in Eqs. (|30|) and (|3Tj) . these 4x4 matrices form a 

sparse matrix of (non-contiguous) 4x4 blocks in the A-body problem. Specially coded 
subroutines are used to efficiently perform these operations. 

Most of the time in VMC or GFMC calculations is spent evaluating wave functions (or 
in GFMC making a propagation step which is equivalent). The pair operators dominate 
this time. The evaluation of a kinetic energy involves numerical second derivatives which 
require 6 A wave function computations. Hence the product of A 7 the number of pairs, and 
of the length of the spin-isospin vector is a good indication of how the total computational 
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Scaling of 


wave function computation 


time 



time scales with A. Table U shows this scaling for various nuclei, assuming M=0 for even 
J nuclei and that good-isospin bases arc being used. The final column shows the product 
of the first three columns relative to 8 Be. We can do calculations up to A=10 routinely 
and a few 12 C calculations have been done. It is clear that this approach is not reasonable 
for 16 O. The last two lines are for "neutron drops" for which isospin does not have to be 
considered. This allows somewhat bigger A to be reached. 

4'3. A Variational Monte Carlo Calculation. - The basic steps in a variational calcu- 
lation are 

• Generate a random position: R = f\, r%, ■ ■ ■ , Ta ■ 

• Make many (1000's) random steps based on the probability P = ^^(R)! 2 . 

• Start integration loop: 

- Make order 10 steps based on P . 

- Compute and sum i? local (R) = [* T (R) t i/* T (R)]/|* T (R)| 2 . Gradients and 
Laplacians are computed by differences: QA evaluations of 1 J , t(R + ^j^) ■ 

• (* T |i/|* T )/(* T |* T ) = average( J ff local ) . 

A random step from a given position, R, to a new position R' is made using the 
Metropolis method: 
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Fig. 5. - Comparison of VMC and GFMC energies for the AV18+IL2 Hamiltonian. The light 
shading shows Monte Carlo statistical errors. 



• Use 3A uniform random numbers on (0,1), {vjj}, to make AR; Axj = 2Sr(uij — 1) 

• Set R' = R + AR, and compute P(AR) = |* t (R')| 2 /I*t(R)| 2 • 

• Make another random number on (0,1): p 

• If P > p, the step is accepted; replace R with R' . 

if P < p, the step is rejected; discard R' and stay at R . 



4'4. Accuracy of VMC energies. - Figure [5] compares VMC energies of various nuclear 
states with the corresponding GFMC values for the AV18+IL2 Hamiltonian. As is de- 
scribed in the next section, the GFMC results are believed to be accurate to 1-2%. For 
4 He the VMC result is quite close to the GFMC. However as we move into the P shell, 
the VMC results get steadily worse. In fact, although the GFMC calculations show that 
this Hamiltonian binds the nuclei shown with the exception of 8 Be, the VMC energies 
are all above the VMC energies for the subclusters that the nuclei can breakup into. 
Furthermore the 8 Be VMC energy is actually lower than those of 9- and 10-body nuclei. 
Calculations with simpler Hamiltonians show that these failures of the VMC energies are 
related to the tensor force; VMC calculations for the simple AV4' potential discussed in 
Sec. 12 31 are quite accurate, while those for AV6' have significant errors. 
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5. Green's Function Monte Carlo General Description 

As shown above, our VMC trial wave functions are not good enough for P-shell nuclei. 
This means that they contain admixtures of excited-state components in addition to the 
desired exact ground-state component, 

(34) * T = * + ^a^i . 

We use Green's Function Monte Carlo to project out of by propagating in imag- 
inary time, t: 

(35) *(t) = cxp{-(H - E )t]*t , 

(36) = e -( E °-^ T x [* + J2 a i e ~ {Ei ~ E ° )T *i\ - 

(37) lim *(r) oc * , 

r— »oo 

where E is a guess for the exact Eq. 

The eigenvalue Eq is calculated exactly while other expectation values are generally 
calculated neglecting terms of order — ^t\ 2 and higher. In contrast, the error in the 
variational energy, Et, is of order — ^t\ 2 , and other expectation values calculated 
with have errors of order | ^ — | • 

The evaluation of *(r) is made by introducing a small time step, At, t = nAr, 



(38) *(r) 



e -(H-E )^ 



where G is the short-time Green's function. The 'J'(t) is represented by a vector function 
of R, and the Green's function, G Q( g(R',R) is a matrix function of R' and R in spin- 
isospin space, defined as 

(39) G a/3 (R',R) = (R',a| e - (ff - Bo)AT |R,/3) . 

It is calculated with leading errors of order (At) 3 as discussed below. Omitting spin- 
isospin indices for brevity, *(R„,t) is given by 

(40) *(R„, t) = J G(R„, R n _i) • • • G(Ri, Ro)tfr(Ro) , 
and 

/ <4(R„) Gt(R n ,R n _ 1 )...Gt(R 1 ,R ) H * T (R ) dP 



(41) E(t) 



/^(R„) Gt(R„,R„_ 1 )..-Gt(R 1 ,R ) * T (R ) dP 



where dP = c?R dRi • • • rfR„. Here we have placed the ^(t) to the left side of H because 
the derivatives in H may be evaluated only on ^t] we cannot compute gradients or 
Laplacians of ^(t). This 3yln-dimensional integral is computed by Monte Carlo. 
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5T. The Short-Time Propagator. - The success of a GFMC calculation depends on 
an accurate and fast evaluation of the short-time propagator, G Q| g(R', R). One wants to 
be able to do this for the largest possible value of At to reduce the number of steps, n, 
needed to reach some asymptotic value of r. The most important features of ^(t) are 
induced by the NN potential, so consider first G a p(R' ,R) for a Hamiltonian with no 
NNN potential. This can be written as 



(42) 



G Q/3 (R',R) = e BoAT G (R',R)(a 



1/3} , 



where 

(43) G (R, R) = (R'|e~ KAr |R) 



2nh 2 AT 



3,1 



exp 



-(R-R) 5 



2h 2 Ar/m 



is the many-nucleon free propagator and go^ij is the corresponding two-nucleon free prop- 
agator, 



(44) 



1 3 



2tt/i2At 



exp 



2n 2 Ar//i 



and \i = m/2 is the reduced mass. 

The Go(R',R) is included in the Monte Carlo integration [Eq. (|4"Tj) ] by using it to 
make the step from R to R'. The magnitudes of the 3 A steps (x, y, and z for each 
nucleon) are determined by sampling a Gaussian of the width given in Eq. (|43[) and the 
directions of the steps are picked by importance sampling; see Ref. |13) for details. 

Eq. (|42f introduces the exact two-body propagator, 



(45) 
(46) 



9ij( r ij > r ij) — { r ij\ e !J \ r ij) J 
ft 2 

Hij = Vfj + Vij . 



All terms containing any number of the same z>ij and K are treated exactly in this 
propagator, as we have included the imaginary-time equivalent of the full two-body 
scattering amplitude. Eq. (|4"2"| still has errors of order (At) 3 , however they are from 
commutators of terms like VijKvik(Ar) 3 which become large only when both pairs ij 
and ik are close. 

To calculate <?y, we use the techniques developed by Schmidt and Lee [55] for scalar 
interactions. These allow <?y to be calculated with high (~ 10 digit) accuracy. However, 
this calculation is quite time consuming. Therefore, prior to the GFMC calculation, we 
compute and store the the propagator on a grid. For a spin-independent interaction, 
the propagator g.ij would depend only upon the two magnitudes r' and r and the angle 
cos(8) = f ' • f between them. Here, though, there is also a dependence upon the spin 
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quantization axis. Rotational symmetry allows one to calculate the spin-isospin compo- 
nents of gij (r',r) for any r' and r by simple SU3 spin rotations and values of gij on a 
grid of initial points r = (0,0, z) and final points r' = (x',0, z'). In addition, the fact 
that the propagator is Hcrmitian allows us to store only the values for z > z' . 

Returning to the full Hamiltonian including NNN forces, the complete propagator is 
given by 

G a/3 (R',R) = e ^G (R,R)exp[-]T(V^(R) + Vf fe (R))^] 



(47) (a|7 3 (R')l7)(7l 
with 

(48) I 8 (R) = 



g j~J 9ij ( r ij ' r ij ) 



2 

\6)(6\I 3 (K)\P) 



The exponential of V?%. is expanded to first order in At thus, there are additional error 
terms of the form V^V^J, k , (At) 2 . However, they have negligible effect since has a 
magnitude of only a few MeV. It was verified that the results for 4 He do not show any 
change, outside of statistical errors, when At is decreased from 0.5 GeV -1 . 

5'2. Problems with Nuclear GFMC. - While GFMC is in principal exact for the (if), 
there are several practical difficulties that make it only approximate. We have made 
many tests of the accuracy of the GFMC energies, both by comparison to other methods 
and by comparing calculations with different At, starting and other computational 
parameters. These tests show that our results for energy are good to ~1% up to ^2% 
for larger A or N — Z ( 8 He is particularly difficult). Some of the problems are: 

5'2.1. Limitation on if. The exact propagator of Eq. (|4T))) can be computed for the 
full Wig potential, however the L 2 and (L • S) 2 terms in the potential correspond to state- 
dependent changes of the mass appearing in the free Green's function. Since we do not 
know how to sample such a free Green's function, we cannot use the exact for the full 
potential, but rather must use one constructed for an approximately equivalent potential 
that does not contain quadratic L terms, namely the AV8' introduced in Sec. 12 31 The 
difference between the desired and approximate potentials is computed perturbatively. 
Comparisons with more accurate calculations for 3 H and 4 He suggest that this introduces 
errors of less than 1%. 

5'2.2. Fermion sign problem. The G(Ri,Rj_i) is a local operator and can mix in 
the boson solution. This has a (much) lower energy than the fermion solution and thus 
is exponentially amplified in subsequent propagations. In the final integration with the 
antisymmetric $t, the desired fcrmionic part is projected out in Eq. (|41[) . but in the 
presence of large statistical errors that grow exponentially with t. Because the number 
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of pairs that can be exchanged grows with A, the sign problem also grows exponentially 
with increasing A. For A>8, the errors grow so fast that convergence in r cannot be 
achieved. 

For simple scalar wave functions, the fermion sign problem can be controlled by not 
allowing the propagation to move across a node of the wave function. Such "fixed- 
node" GFMC provides an approximate solution which is the best possible variational 
wave function with the same nodal structure as ^t- However, a more complicated 
solution is necessary for the spin- and isospin-dependent wave functions of nuclei. This 
is provided by "constrained-path" propagation in which those configurations that, in 
future generations, will contribute only noise to expectation values are discarded. If the 
exact ground state \^o) were known, any configuration at time step n for which 

(49) *(R„)t* (R„) = , 



where a sum over spin-isospin states is implied, could be discarded. The sum of these 
discarded configurations can be written as a state \^d), which obviously has zero overlap 
with the ground state. The "i>d contains only excited states and should decay away as 
r — > oo, thus discarding it is justified. Of course the exact is not known, and so 
configurations are discarded with a probability such that the average overlap with the 
trial wave function, 

(50) (* d |* r ) = 0. 

Many tests of this procedure have been made and it usually gives results that 
are consistent with unconstrained propagation, within statistical errors. However a few 
cases in which the constrained propagation converges to the wrong energy (either above 
or below the correct energy) have been found. Therefore a small number, n u = 10 to 
20, of unconstrained steps are made before evaluating expectation values. These few 
unconstrained steps, out of typically 400 total steps, appear to be enough to damp 
out errors introduced by the constraint, but do not greatly increase the statistical error. 
Unfortunately, the constrained-path E(t) are not upper bounds to the true Eg; examples 
have been found in which the constrained energies evaluated with inadequate n u are below 
E . 

5'2.3. Mixed estimates extrapolation. As shown in Eq. (|4"Tj) . GFMC computes 
"mixed" expectation values between \&y and ^(r) of operators, 

(51) ( ° >Mlxcd - <*(r)|* T ) ■ 

The desired expectation values, of course, have ^(r) on both sides. By writing \I/(t) = 
^T + S^ir) and neglecting terms of order [(^(r)] 2 , we obtain the approximate expression 

( 52 ) (0(r)) = \^)\*(t)) « (0(r))Mixcd + [(0(r)) Mixcd - (Q) T ] , 
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(0) T 


<°>Mix 


(0} Mix " (Oh 


<o> 


K 


146.60 


153.49 


6.89 


160.39 


vmic 


-171.64 


-180.49 


-8.85 


-189.34 


vc 


1.53 


1.54 


0.01 


1.55 


Vijk 


-4.10 


-6.43 


-2.34 


-8.77 


Sum 


-27.61 


-31.90 


-4.28 


-36.17 


H 


-27.61 


-31.90 




-31.90 



Table II. - Contributions to {H} for 6 Li (MeV) 



where (0}t is the variational expectation value. More accurate evaluations of (0(r)) are 
possible, essentially by measuring the observable at the mid-point of the path. However, 
such estimates require a propagation twice as long as the mixed estimate and require 
separate propagations for every (O) to be evaluated. 

The expectation value of the Hamiltonian is a special case. The (H (r))Mixod can be 
re-expressed as [23] 

^ T \ e -i,H-E 0) r/2 He - { H-E a) r,2^ T) 
(53) ( J ff(r)) Mi xcd - ^ r | e _( H _ Bo ) T/2e -(if-£ ) r /2|^ r ) ^ E » > 

since the propagator exp[— (H— Eq)t] commutes with the Hamiltonian. Thus (-ff(r)) Mixed 
is already the correct expectation value and must not be extrapolated. This results in 
the unfortunate circumstance that the sum of the pieces of (H) is not equal to the full 
GFMC value of (H). An example is shown in Table HT1 which shows VMC, mixed, and 
extrapolated energies for 6 Li computed with the AV18+IL2 Hamiltonian. The sum of 
the extrapolated kinetic and potential energy values is 4.3 MeV different from the total 
energy. This means that the extrapolated values of the pieces have errors whose absolute 
sum is at least this big. 

Instead of the linear extrapolation of Eq. (|52[) . one can also use a ratio extrapolation: 



(54) (0(r)> « fefeed 

{U)t 

which is the same as Eq. (|52p to lowest order in [(0(r))Mixcd — (O)t], but obviously has 
different quadratic errors. This method has the feature that if both (0(r))Mixed and 
{0)t have the same sign, then the extrapolated (0(t)) will also have that sign. This 
is an advantage for quantities such as densities (see Sec. [8]) which must be positive; if 
the GFMC is reducing the density at large r where it is exponentially falling, linear 
extrapolation can result in negative values. Of course such large extrapolations by either 
method are uncertain. 
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5'3. A Simplified GFMC Calculation. - The basic steps in a GFMC calculation are 

• Start with collection of 4 r (r = 0, Rj) = ^t(Rj) from a VMC calculation 

• Loop over time steps t„ = rtAr 

— Loop over configurations j 

* Make a random step to R'- = Rj + ARj by sampling Gq(R', R) 

* Sample several directions based on simplified and potential 

* Compute *(r n+ i,Rj) = G(Rj, Kj)^(r n , Rj) 

* Possibly mark as killed due to the constraint Wfrn+i, R'j) ■ ^T(Rj) 

* Use importance sampling to kill or replicate the configuration ty(T n +i, R^) 

— Every 20-40 time steps 

* Compute the local energy ^(r n , Kj)H^! T (ILj)/^(T n , Rj)\& T (Rj) 

* Check that total number of configurations is staying reasonably constant 

GFMC calculations are quite computer intensive. For example, a typical 8 Li calcu- 
lation requires 300 processor hours running at a delivered (not theoretical-peak) speed 
of one GigaFLOPS. As shown in Table HI a 10 B calculation will need about 10 times 
this and 12 C 250 times it. Clearly such calculations are practical only on highly parallel 
computers. Our GFMC program uses a master-slave structure in which each slave gets 
a number of configurations to propagate as outlined in the preceding paragraph. The 
computed energy results are sent back to master for averaging as they are generated. 

Because configurations are multiplied or killed during propagation, the work load 
fluctuates. It is important to periodically rebalance the work load - otherwise slaves will 
wind up with nothing to do while the last slave with the biggest work load finishes its 
calculations. To do this, the master periodically collects load statistics and then tells 
slaves to redistribute some of their configurations. The slaves have work (energy calcu- 
lations left from previous time steps) set aside to do during this synchronization. This 
method results in parallclization efficiencies of typically 95% on up to 2000 processors. 
However, the next generation of large computers will have 10,000 to 100,000 processors, 
which is more than the number of configurations to be propagated for a large nucleus 
like 12 C. Thus the program has to be made parallel at a finer level; this is being worked 
on. 

5'4. Examples of GFMC propagation. - Figure [5] shows the E{t) as a function of the 
imaginary time, t, for the beginning of GFMC propagation for i He. The propagation 
starts from the VMC value of -26.92 MeV, and initially decreases rapidly with increasing 
r; essentially converged values are achieved by r = 0.01 MeV -1 . The propagation is 
continued another factor of 10 to r = 0.1 MeV -1 and the results averaged over the last 
half of the propagation to get the converged result of -28.300(15) MeV. The curve is a 
fit of the form 

(55) E(r) = E > ^Ete^-Etr) 



l + E^ezK-^r) 




to the computed E(r) using three terms. The E* was fixed at the first + excitation 
energy of 20.2 MeV of 4 He, and the other two E* and the three oti were varied in the 
fit. The fitted E* turn out to be very large, 340 and 1480 MeV, with small a u 0.0018 
and 0.00046, respectively. Thus the errors in the VMC ^>t correspond to small amounts 
of extremely high excitation energy; GFMC is particularly efficient at filtering out such 
errors. 

Figure [7] shows GFMC propagation, using the AV18+IL2 Hamiltonian, of the ground, 
first 3 + , and 2 + states of 6 Li. The propagation for the ground state (which is particle 
stable with this H) and the 3 + (which is only slightly above the d+a threshold and 
experimentally has a narrow width) is stable after r = 0.2 MeV -1 . However the 2 + (a 
broad resonance) never becomes stable; the E(t) are decaying to the threshold energy of 
separated a and d clusters. Because the 3 + state E(t) stops decreasing around r=0.2, 
the £'(r=0.2), shown by the star, is best GFMC estimate we can currently make of 
the resonance energy. However it is now possible to make GFMC calculations using 
scattering-wave boundary conditions (see Sec. [7|) and this method will be applied to 
states such as 6 Li(2 + ). 



6. Results for energies of nuclear states 

Figure [3] compares the GFMC energies of various nuclear states with experiment and 
Fig. [S] does the same for excitation energies. In both cases the left set of bars for each 
isotope shows results using just the AV18 NN potential while the middle set of bars is for 
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Fig. 7. - GFMC propagation for three states of 6 Li. 



2.5 



the full AV18+IL2 Hamiltonian. As has already been observed, AV18 alone significantly 
underbinds all nuclei except the deuteron; including IL2 results in fairly good agreement 
with the experimental values. The excitation spectra in Fig. [8] show that IL2 also fixes 
other problems that arise when just a NN potential is used. For example, spin-orbit 
splittings are usually too small without the WWW potential (note the | — | and 
| — | splittings in 7 Li and the \ — | splitting in 9 Li). As is discussed in the next 
subsection, even the ordering of states can be changed by the NNN potential. 

The discussion of Sec. [5] implies that GFMC can be used only for the lowest state 
of each set of quantum numbers but Fig. [3] shows several states with the same J*'. The 
ability of GFMC to provide such results was demonstrated in Ref . [ST] . 

6T. Ordering of States in 10 Be and 10 B. - Figure [9] shows the beginning of the com- 
puted and experimental excitation spectra of 10 Be and 10 B. We see that NN potentials 
with no NNN predict a 1 + ground state for 10 B while the Illinois-2 NNN potential fixes 
this and gives the correct 3 + ground state. No-core shell model calculations show that 
other AW potentials without NNN potentials also give a 1 + ground state for 10 B |24) . 
so this is not a failure of just AV18. This incorrect ground-state prediction is another 
manifestation of too-small spin-orbit splitting using just WW potentials; in 1956 D. Ku- 
rath showed that the relative positions of the 3 + and 1 + levels depends on the amount 
of spin-orbit strength in a shell-model calculation 25j. 
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Fig. 8. - GFMC computations of excitation energies for the AV18 and AV18+IL2 Hamiltonians 
compared with experiment. The shading on the experimental energies shows the widths of 
resonances. 



The first two excited states in 10 Be are both 2+ and Fig. [H] shows that including the 
IL2 NNN potential reverses their order. VMC and GFMC calculations predict large 
positive and negative quadrupole moments (Q) for these states; with no NNN potential, 
the GFMC energy of the Q > state is the lower, while adding the IL2 NNN potential 
reverses this. VMC also predicts a large B(E2) to the g.s. for only the state with Q > 0. 
An ATLAS experiment for the B(E2) and quadrupole moments of these states will be 
made to determine if the reversal of order given by IL2 is correct. 

6'2. Charge Dependence and Isospin Mixing. - The differences of the energies of states 
in the same isomultiplet are a probe of isospin-breaking components of the Hamiltonian. 
The largest such component is the Coulomb potential between protons, but it has been 
known since 1969 that this does not fully account for the measured differences [55]; the 
discrepancy is referred to as the Nolen-Schiffer anomaly. It is convenient to parametrize 
the energies of the isomultiplets by the coefficients a^ T , 

(56) E A , T (T Z ) = a ( ^ T Q n {T,T z ) , 

n<2T 
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where Qo = 1 is the isoscalar component, Q\ — T z the isovector, and Q2 = |(3T| — T 2 ) 
the isotensor. Term 18 of AV18 contributes to the isovector component and terms 15 
to 17 to the isotensor component; both terms receive contributions from the various 
electromagnetic terms in AV18. Table ITTll shows computed and experimental values of 
the coefficients for several isomultiplets; the ^ CSB + CD column gives contribution of the 
nuclear (strong-interaction) terms in AV18, while the K column shows that resulting 
from the difference of the proton and neutron masses. In general the non-Coulomb 
electromagnetic and strong charge symmetry breaking and charge dependent terms result 
in good agreement with the experimental values and thus resolve the Nolen-Schiffer 
anomaly. 

Figure ITD1 shows the beginning of the experimental excitation spectrum of 8 Be. From 
16 to 19 MeV excitation there are pairs of 2 + , 1 + , and 3 + levels with isospin and 1. 
The two 2 + levels are very close and hence strongly isospin mixed; the 3 + levels also 
have significant mixing. The mixing has been experimentally known from the decay 
properties of the states since 1966 (27]. However, as with the Nolen Schiffer anomaly, 
calculations using just Coulomb mixing underestimate the amount of mixing. Table IIVI 
shows GFMC calculations using AV18+IL2 of the mixing matrix elements for the three 
pairs [28) . The contribution of the nuclear CSB term is relatively more important here 
than for the Nolen Schiffer anomaly. The agreement of the predictions with the data is 
not as satisfactory as for the Nolen Schiffer anomaly. The final line of the table shows the 
mixing matrix element between the first 2 + (at 3 MeV) state and the isospin-1 17-MeV 
state. The small value of this matrix element, and the large energy denominator, shows 
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T 


n 


^Coul 


otherEM 


^CSB+CD 


^-CSB 


Total 


Expt. 


3 H- 3 He 


i 

2 


1 


649 


29 


64 


14 


757 


764 


7 Li- 7 Be 


1 

2 


1 


1458 


40 


83 


23 


1605 


1644 


7 He, 7 Li*, 7 Be*, 7 B 


3 
2 
3 
2 


1 

2 


1286 
132 


14 
7 


49 
34 


17 


1366 
174 


1373 
175 


8 Li, 8 Be*, 8 B 


1 

1 


1 

2 


1692 
140 


24 
5 


78 
-5 


24 


1818 
140 


1770 
145 


8 He, 8 Li*, 8 Be*, 8 B*, 8 C 


2 
2 


1 
2 


1719 
153 


13 
7 


83 
42 


26 


1840 
203 


1659 
153 



Table III. - Computed (AV18+IL2) and experimental Nolen-Schiffer energies for several iso- 
multiplets 



that the first 2 + state has very little T = 1 contamination; this is important to the 
possibility of using 8 Li(/3 - ) 8 Be(l st 2 + ) decay as a test of V — A. 

6'3. Can Modern Nuclear Hamiltonians Tolerate a Bound Tetraneutron? . - In 2002 a 
claimed observation of a bound tetraneutron was published [211 [3D] ■ The experiment did 
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Fig. 10. - Experimental spectrum of Be 
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J p GFMC Expt 

Coulomb Strong CSB Other Total 

2 nd 2+ 78 21 16 115 144 

1+ 80 18 4 102 120 

3+ 61 15 14 90 63 

1 st 2+ 4 0.4 1 6 



Table IV. - Isospin mixing matrix elements for Be in keV 



not produce a definite binding energy, just the statement that four neutrons are weakly 
bound. It is well known that the dineutron is not bound, but rather has low-energy pole 
on the second sheet (pseudo bound state); AV18 along with other realistic potentials has 
this feature. A set of GFMC calculations were made to see if the AV18+IL2 Hamiltonian, 
or acceptable modifications of it, could reproduce the tetraneutron claim [31] ■ It was 
clearly established that the unmodified AV18+IL2 does not bind 4 n; at most there is 
some weak resonance at E ~ +2. MeV. 

Minimal modifications to AV18+IL2 to give E( 4 n) ~ —0.5 MeV were then made 
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Fig. 11. - Nuclear energies with Hamiltonians modified to bind a tetraneutron 
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Fig. 12. - One-body (left) and two-body (right) densities of tetraneutrons. Circles are for the 
AW and squares for the T = § NNN modifications; the diamonds are for 2 n. 



and the effects of such modifications on the energies of other, well established, nuclei 
were computed. Figure [TT1 shows some of these. In the first case (the left-hand bars), 
the intermediate-range part of the ^Sq partial-wave potential in of AV18 was increased 
enough to bind four neutrons. This increase results in the dineutron also being bound, 
in fact the 4 n is still unbound against breakup into two dineutrons in this model. As the 
figure shows, other existing nuclei ( 3 H, 4 He, etc.) all become significantly over bound 
with this increased potential. Also 6 n and 8 n (not shown) become bound. An attempt 
to bind the 4 n by changing the 3 Pj part of AV18 was also made, but this requires a huge 
change which very strongly over binds other nuclei. 

A second attempt was to add an attractive Vijk(T = |) to H. This has two advan- 
tages: 1) it has no effect on NN scattering and does not make a bound 2 n; 2) because 
the modification is made only in isospin-| triples, it has no effect on 3 H, 3 He, or 4 He; 
these have only isospin-^ triples. However, as can be seen in the figure, as soon as this 
potential can act (i.e. in nuclei with T = | triples), it produces huge overbinding. The 
most dramatic effect is in pure neutron systems: 6 n is bound by 220 MeV and 8 n by 650 
MeV. These are the most stable 6- and 8-nucleon systems with this Hamiltonian, so all 
other 6- and 8-body nuclei would beta decay to them! 

Finally an attractive four-nucleon, T — 2, potential was added (not shown). This 
does not effect 6 Li but does very strongly overbind 6 He, larger nuclei, and pure neutron 
systems with more than four neutrons. The conclusion of the study is that a bound 4 n is 
incompatible with our understanding of nuclear forces. In the meantime the experiment 
has not been successfully reproduced. 

The % NN and T = § NNN modifications have very different effects on A > 4 
binding energies, even though both have been adjusted to bind 4 n by only 0.5 MeV; 
the NNN modification results in much more severe overbinding. Figure [12] helps to 
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explain this. The left panel shows the one-body densities of the two bound 4 n systems; 
the NN modification results in a very diffuse 4 n with a rms radius of 8.9 fm. The 
right panel shows the pair density which is proportional to finding two neutrons a given 
distance apart. For the 1 Sq NN modification, the density for a S=0 pair is peaked around 
1 fm and is very similar to the pair distribution of the isolated 2 n which is bound by 
this potential. The 5=1 pair distribution is peaked at 12 fm; it arises from neutrons 
in different dineutron clusters. Thus this 4 n looks like two widely separated dineutrons. 
The binding comes from the small tails where neutrons from each dineutron get close 
enough to interact. The change of i>jviv( 1 <S'o) needed to achieve this is not that big and 
hence bigger nuclei overbound only somewhat. 

On the other hand, a Vnnn requires all three neutrons in a triple to be close together 
to be effective. Thus the T = | NNN modification must bring the two dineutrons close 
together. This results in the much more compact one-body density, with a rms radius of 
only 1.9 fm, shown in the left panel and a two-body density that is much sharper than 
the isolated 2 n. The Vnnn(T= §) must be large to achieve this high density and it thus 
has a large effect in all bigger nuclei. 

7. GFMC for Scattering States 

The GFMC calculations presented so far have treated the nucleus as a particle-stable 
system; that is the starting wave function, vf^, is exponentially decaying as any nucleon 
is removed to a large distance from the center of mass. However many of the states of 
interest are particle unstable; they are above the threshold for emission of a single nucleon 
or, as is often the case for light nuclei, the threshold for breakup into subclusters. This 
approximation appears to be adequate for narrow resonances which have only a small 
scattering- wave component. However broad resonances should really be computed with 
proper scattering-wave boundary conditions; as shown in Fig. [7] we do not achieve a 
converged energy for a broad state with the bound-state GFMC. In addition to being 
the correct approach, scattering solutions allow one to compute the phase shift as a 
function of energy and thus obtain the width of the resonance. Finally one might also 
be interested in the phase shifts for partial waves that have no resonance. 

Preliminary GFMC calculations of neutron-alpha scattering were made in 1991 [32], 
but detailed, high statistics results have been obtained only recently [33J . Instead of using 
exponentially decaying wave functions, we construct to have a specified logarithmic 
derivative, 7, at some large boundary radius (R > 7 fm). Here R is the maximum 
distance that any nucleon is allowed to get from the other A — 1 nucleons. The GFMC 
propagation uses a method of images to preserve 7 at R, and thus finds E(R, 7); the 
eigenenergy that corresponds to the boundary condition. The phase shift, 6(E), can then 
be computed from R, 7, and E(R, 7). This procedure is repeated for a number of 7 and 
6(E) is mapped out parametrically. 

Figure fT31 shows an example of this for the broad | resonance in 5 He. The bound- 
state boundary condition docs not give a stable energy, but is decaying to the n+ 4 He 
threshold energy. The scattering boundary condition produces a stable energy; the value 
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Fig. 13. - GFMC propagation for 5 He(^ ) using bound- and scattering-state boundary condi- 
tions (diamonds and circles, respectively). 



of 7 used in this example results in an energy slightly above the resonance energy. 

Figure [T4l compares the calculations of n+a scattering with a i?-matrix analysis of the 
data [33] (solid curves). The left panel shows the partial- wave phase shifts computed for 
three different Hamiltonians: AV18 with no V ijk , AV18+IL2, and AV18+UIX (UIX is an 
older NNN potential [35] that, with AV18, correctly binds 3 H and 4 He, but underbinds 
P-shell nuclei). All three Hamiltonians give very similar results for the i + partial wave 
which has no resonance. However only the AV18+IL2 correctly reproduces the two P- 
wave partial waves; AV18 alone misses both of them and AV18+UIX fits the | partial 
wave but has too-small spin-orbit splitting and misses the | one. The right panel shows 
the partial-wave cross sections for the AV18+IL2 Hamiltonian. The strong resonance in 
the | channel is very well reproduced showing that both the position and width of the 
resonance agree with the data; the much broader i resonance is also well reproduced. 
In addition the good agreement with the low-energy i + cross section data shows that 
the scattering length is reproduced. 

This first study is very promising; the GFMC method, with its ability to have cor- 
rect asymptotic forms, should be applied to other scattering calculations including a 
number of broad resonances [ 7 < 9 He, 6 Li(2+), 8 Be(2+,4+), etc] and the initial states of 
astrophysically interesting capture reactions [ 4 He(d,7) 6 Li, 7 Be(p,7) 8 B, etc]. 
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Fig. 14. - GFMC calculations of n + a scattering in the three principal partial waves. The 
left panel shows partial-wave phase shifts for three Hamiltonians. The right panel shows corre- 
sponding partial-wave cross sections for the AV18+IL2 Hamiltonian. The experimental data is 
represented by the solid curves. 



8. Coordinate- and momentum-space densities 

Up to now we have been concentrating on GFMC calculations of energies of nuclear 
states. However matrix elements of any operator may be evaluated in the GFMC propa- 
gation by using the extrapolation formulas, Eqs. (|o^|) or I discuss some recent work 
on the charge radii, the corresponding densities, and momentum-space densities in this 
section. 

8T. RMS radii and one-body densities of helium isotopes. - A few years ago, a group 
at Argonne measured the RMS charge radius of the radioactive nucleus 6 He (/J-decay 
half-life 0.8 sec.) with the remarkable accuracy of 0.7% [35] and the corresponding 
measurement for 8 He (/3-decay half-life 0.1 sec.) will be published soon [37]. This has 
led us to attempt equally precise GFMC calculations of the corresponding point proton 
RMS radii. Such calculations are very difficult because of the small separation energies 
of the two valence neutrons in these isotopes (i?sep = 0.97 MeV for 6 He and 2.14 MeV 
for 8 He). Changes in the starting and other aspects of the GFMC calculations can 
result in changes of 200 keV (400 kcV for 8 He) in the computed energy (and hence 
Esep)- The RMS radius depends strongly on E'sep! as E'sep goes to zero, the radius 
goes to infinity. Thus we cannot give a precise value for the computed RMS radius for a 
specific Hamiltonian. Instead we find that the computed values for the same Hamiltonian 
with different GFMC calculations, or even for different Hamiltonians, all lie in a band 
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Fig. 15. - GFMC calculations of point proton RMS radii of 6,8 He plotted as a function of the 
two-neutron separation energy obtained in the calculation. 



of radius versus separation energy. This is shown in Fig. IT51 which shows results for two 
Hamiltonians, AV18+IL2 and AV18+IL6, each with several GFMC calculations (IL6 is 
a newer, unpublished, version of IL2). The stars in each panel show the experimental 
point radii at the experimental separation energies; they are clearly consistent with our 
calculations which give 1.92(4) fm for 6 He and 1.82(2) fm, for 8 He. These numbers are 
both significantly bigger than the RMS point radius of 4 He which is 1.46 fm. 

Figure [TBI shows the point proton and neutron densities of 4 ' 6,8 He. The alpha particle 
is extremely compact; its central density is twice that of nuclear matter. In these calcu- 
lations it has identical proton and neutron densities which is a very good approximation. 
As is shown below, the valence neutrons in 6,8 He do not seriously distort the He core, 
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Fig. 17. - GFMC calculations, using AV18+IL2, of two-proton densities of helium isotopes. 



rather they just drag the 4 He center of mass around. This results in the proton density 
being spread out which is why the charge radii of 6,8 He are so much greater than that 
of 4 He even though all three nuclei have just two protons. The right panel of the figures 
clearly shows that 6 ' 8 He have large neutron halos due to the weak binding of the extra 
neutrons. The neutron halo of 6 He is more diffuse than that of 8 He as is expected from 
the smaller -Esep of 6 He. 

8'2. Is an alpha particle in a sea of neutrons still an alpha particle?. - The previous 
subsection showed that the proton density in 6,8 He is much more spread out than the 
density of 4 He, even though 6 ' 8 He have only extra neutrons added to a 4 He core. This 
might first be thought to indicate that the core of 6 ' 8 He has been considerably enlarged 
by the neutrons. This can be studied by computing p pp , the pair density which is propor- 
tional to the probability for finding two protons a given density apart. These distribution 
functions are shown in Fig. [T71 again calculated with GFMC for the AV18+IL2 model. 
These nuclei each have just one pp pair which presumably is in the "alpha core" of 6 ' 8 He. 
Unlike the one-body densities, these distributions are not sensitive to center of mass 
effects, and thus if the alpha core of 6 ' 8 He is not distorted by the surrounding neutrons, 
all three p pp distributions in the figure should be the same. 

We see that the pp distribution spreads out slightly with neutron number in the helium 
isotopes, with an increase of the pair rms radius of approximately 4% in going from He 
to 6 He, and 8% to 8 He. While this could be interpreted as a swelling of the alpha core, it 
might also be due to the charge- exchange (tj • Tj ) correlations which can transfer charge 
from the core to the valence nucleons. Since these correlations are rather long-ranged, 
they can have a significant effect on the pp distribution. VMC calculations of 4 He with 
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Fig. 18. - Two-nucleon momentum distributions in 4 He, computed by VMC for the AV18+UIX 
Hamiltonian. The symbols show pp distributions; the curves are for np. 



wave functions modified to give p pp distributions close to those of 6 ' 8 He suggest that the 
alpha cores of 6,8 He are excited by ~ 80 and ~ 350 keV, respectively, which corresponds 
to only a 0.4—2% admixture of the first + excited state of 4 He at 20 MeV. Thus almost 
all of the increased RMS radius of the proton density is due to the a core of 6 ' 8 He being 
pushed around by the neutrons and not distortions of the core. 

8'3. Two-nucleon knockout - (e,e'pN). - A recent JLAB experiment for 12 C(e,e'pN) 
measured back to back pp and np pairs; that is pairs with total CM. momentum Q = 0, 
as a function of their relative momentum, q |38j . They found that the cross section for np 
pairs with q in the range 2-3 fm -1 is 10-20 times larger than that for pp pairs in the same 
range. To study this we made VMC calculations of the corresponding pair momentum 
distributions (pnn) m several nuclei from 3 H to 8 Be [39]. The calculations for 4 He are 
shown in Fig. [TH1 Results for the AV18+UIX Hamiltonian (which for ^4=3, 4 is a good 
approximation to AV18+IL2) are shown as the solid line (np pairs) and solid circles (pp 
pairs). Around q=2 fm -1 there is a deep minimum in the pp density which results in 
large values for p n p/p P p- The dashed curve and open circles show the corresponding 
quantities computed for the AV4' AW potential (Sec. I2'3[) with no NNN potential; in 
this case both densities have a deep dip and there is no enhancement of the ratio. 

The AV4' potential has no tensor force and thus np pairs are just S'-wave while, for 
the full Hamiltonian, isospin-0 np pairs, like the deuteron, have a D-w&ve admixture. 
The S'-wave deuteron momentum distribution has a zero at 2 fm -1 which is filled in by 
the D-wave contribution, the same as is seen here for np pairs in He. The tensor force 
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is much smaller in the pure T—l pp pairs, so the deep S- wave minimum is not filled in 
for those pairs, even with the AV18+UIX Hamiltonian. Calculations for 3 He, 6 Li, and 
8 Be all show this effect although the deep minimum is somewhat filled in for 8 Be. Thus 
the JLAB experiment shows the importance of tensor correlations up to > 3 fm -1 . 

9. Conclusions 

Quantum Monte Carlo methods are powerful tools for studying light nuclei with 
realistic nuclear interactions. Calculations of A — 6 - 12 nuclear energies with accuracies 
of 1 — 2% are possible and the AV18+IL2 reproduces binding energies with an average 
error of order 0.7 MeV for A = 3 — 12. The NNN potential is required for overall P-shell 
energies and for spin-orbit splittings and several level orderings. 

The QMC methods allow matrix elements of many operators of interest to be com- 
puted. This contribution presents rms radii, one- and two-body densities and two-body 
momentum distributions. These are generally in good agreement with experiment. Re- 
cently GFMC values of A=6,7 electromagnetic and weak transitions have also been com- 
puted; these improve on older VMC calculations and also generally agree with exper- 
iment [40 . Another topic not covered here is overlap functions and the related spec- 
troscopic factors; these are used as input to calculations (such as distorted-wave Born 
approximation) of nuclear reactions; recent results are presented in Rcfs. (4TJ [42]. A 
just-finished interesting study used VMC calculations to investigate the effects on nuclear 
binding energies of changes in the Hamiltonian induced by changes of the fundamental 
constants [33]. 

GFMC calculations of are very computer intensive and at present 12 C can just barely 
be done. However a new generation of extremely parallel computers is becoming available 
and we are working with computer scientists to enable the GFMC program to make use 
of these machines. This should lead to the possibility of detailed studies of 12 C including 
second + (Hoyle) state which is the doorway for triple-alpha burning. This state has 
resisted precise calculation by shell-model based methods; we hope that our more flexi- 
ble variational wave functions, combined with GFMC propagation, will overcome these 
difficulties. 

But perhaps the most important advance in nuclear GFMC is the computation of 
scattering states. In these calculations the correct scattering-wave boundary condition is 
achieved. The resulting wave functions will be used to compute reactions of astrophysical 
interest such as 3 He+a -^ 7 Be, p+ 7 Be — > 8 B, and n+(a+a) -+ 9 Be. Indeed all big-bang 
nucleosynthesis, solar neutrino, and some r-process seeding reactions are accessible. 

* * * 

As can be seen from the author lists of the citations, the work reported in this con- 
tribution is the result of long-term collaborations with Joesph Carlson (who invented 
nuclear GFMC), Kenneth M. Nollett (who is doing the GFMC scattering), Vijay R. 
Pandharipande (who for many years guided and inspired our group), Rocco Schiavilla 
(who has provided the expertise on the electroweak currents) and Robert B. Wiringa 
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(who has developed the VMC wave functions). I thank R. B. Wiringa for also making 
a critical reading of this MS. The work would not have been possible without extensive 
computer resources provided over the years by Argonne's Mathematics and Computer 
Science Division (most recently on the IBM Blue Gene), Argonne's Laboratory Comput- 
ing Resource Center, and the U.S. Department of Energy's National Energy Research 
Scientific Computing Center. This work is supported by the U.S. Department of Energy, 
Office of Nuclear Physics, under contract DE-AC02-06CH11357. 
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